Thermodynamics of ultrasmall metallic grains in the presence of pairing and exchange 

correlations: mesoscopic fluctuations 
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We study the mesoscopic fluctuations of thermodynamic observables in a nanosized metallic grain 
in which the single-particle dynamics are chaotic and the dimensionless Thouless conductance is 
large. Such a grain is modeled by the universal Hamiltonian describing the competition between 
exchange and pairing correlations. The exchange term is taken into account exactly by a spin- 
projection method, and the pairing term is treated in the static-path approximation together with 
small-amplitude quantal fluctuations around each static fluctuation of the pairing field. Odd-even 
particle-number effects induced by pairing correlations are included using a number-parity projec- 
tion. We find that the exchange interaction shifts the number-parity effects in the heat capacity 
and spin susceptibility to lower temperatures. In the regime where the pairing gap is similar to or 
smaller than the single-particle mean level spacing, these number-parity effects are suppressed by 
exchange correlations, and the fluctuations of the spin susceptibility may be particularly large. How- 
ever, for larger values of the pairing gap, the number-parity effects may be enhanced by exchange 
correlations. 

PACS numbers: 74.78.Na, 74.25.Bt, 75.75.-c, 05.30.Fk 



I. INTRODUCTION 

The physics of nano-sized metallic grains has attracted 
much attention following a series of experiments by 
Ralph, Black and Tinkham^r— in which individual en- 
ergy levels of ultra-small aluminum grains were resolved 
by single-electron-tunneling spectroscopy. Recent ad- 
vances have made it possible to achieve better control 
over the size of the grain, which is important for quanti- 
tative measurements. In the experiments of Ref. very 
high-quality spectra of chemically synthesized gold nano- 
particles were obtained. 

Typical grains used in the spectroscopic experiments 
are in the ballistic regime, i.e., their size is smaller than 
the mean free path, and electron transport is determined 
by scattering from the boundaries of the grain rather 
than from impurities. When the boundaries are suffi- 
ciently irregular, the single-particle dynamics are chaotic. 
This induces sample-specific fluctuations of observables, 
and the meaningful quantities are the statistical distri- 
butions of these observables; see Ref. d and references 
therein. The single-particle energies and wave functions 
follow the statistics of the random-matrix theory (RMT) 6 
in a Thouless energy window i?Th around the Fermi en- 
ergy, where -&rh is determined by the time it takes for 
an electron to move across the grain. 

When .&rh is much larger than the single-particle mean 
level spacing 8, the grain is described by the so-called uni- 
versal Hamiltonian. 7,8 This Hamiltonian contains three 
interaction terms: a "classical" charging energy term, a 
pairing term that is characterized by a bulk pairing gap 
A, and exchange term that depends on the total spin of 
the grain and is characterized by a coupling constant J s . 
These three interaction terms are universal, i.e., they are 
independent of the particular realization of the single- 



particle Hamiltonian. Here we assume J s /8 < 1 so that 
we are below Stoner instability of macroscopic polariza- 
tion. 

When the pairing term is suppressed, (i.e., when 
only charging and exchange terms contribute), thermo- 
dynamic observables of the universal Hamiltonian can 
be calculated in closed form using a spin-projection 
method^ In Refs. [13 and [TJ, a Hubbard-Stratonovich 
transformatio n 12 ' 13 was employed to calculate in closed 
form observables such as the tunneling density of states 
and spin susceptibility. 

In the absence of the exchange term, the univer- 
sal Hamiltonian has the form of the Bardeen-Cooper- 
Schrieffer (BCS)^ Hamiltonian. In the bulk limit A/6 > 
1, an attractive pairing interaction leads to supercon- 
ductivity. Effects of the BCS interaction in nano-sized 
metallic grains were studied extensively; see Ref. [H and 
references therein. Anderson argued^ that the smallest 
possible size of a system that can be a superconductor 
is determined by the condition A/ 6 ~ 1. In the exper- 
iments of Ref. l|4j, a pairing gap was clearly observed 
in the excitation spectra of the largest aluminum grains 
containing an even number of electrons, while it was 
impossible to resolve such a gap in the smaller grains. 
This, however, does not necessarily mean that pairing 
correlations disappear in the smaller grains. It was pro- 
posed that thermodynamic properties could be a more 
suitable tool to probe this fluctuation-dominated regime, 
in which A/ 8 < Signatures of pairing correlations 
in this regime are the dependence of observables on the 
number parity of electrons in the grain. A good exam- 
ple is the re-entrant behavior (i.e., a local minimum) of 
the spin susceptibility with decreasing temperature in an 
odd graini 17 ! 18 Odd-even effects in the heat capacity and 
magnetic susceptibility were experimentally observed in 



2 



small palladium clusters J£ 

BCS theory breaks down when A/S < 1 and fluctua- 
tions of the gap order parameter (beyond its mean-field 
BCS value) are important. In the static-path approxima- 
tion (SPA)j22r— only static fluctuations of the gap are 
taken into account. A better approximation, the SPA 
plus random-phase approximation (RPA), takes into ac- 
count small-amplitude time-dependent quantal fluctua- 
tions of the order parameter around each static fieldi22r— 
Number-parity effects can be studied by using an ex- 
act number-parity projection.— ~— The heat capacity and 
spin susceptibility of a metallic grain (without exchange 
correlations) as functions of temperature were studied in 
the SPA+RPA method together with a number-parity 
projection in Ref. HU as well as by quantum Monte Carlo 
method a 33 i 34 and by Richardson's solution! 17 ' 34 ~ — In all 
of those calculations, number-parity effects were clearly 
identified in both the heat capacity and spin susceptibil- 
ity of the grain. Signatures of pairing correlations were 
also found in the spin susceptibility as a function of mag- 
netic fteldM^S. 

The exchange interaction competes with the BCS-likc 
pairing interaction. Exchange tends to maximize spin 
polarization, while pairing correlations tend to minimize 
the spin. It is known that, depending on the values of 
A/ 5 and J s /S, the ground state of a system can be su- 
perconducting, ferromagnetic, or one in which pairing 
and ferromagnetic correlations coexisti 39 ' 40 The effects of 
mesoscopic fluctuations on this competition were studied 
in Ref. El 

The effect of both pairing and exchange correlations 
on the thermodynamic properties of the grain (heat ca- 
pacity and spin susceptibility) was studied in Ref. |42| for 
the case of an equally spaced single-particle spectrum by 
using a quantum Monte Carlo method. These thermody- 
namic quantities can also be calculated directly from the 
eigenvalues of the universal Hamiltonian using Richard- 
son's solution, modified to take into account the exchange 
interaction^ The combined effect of exchange and pair- 
ing interactions on the spin susceptibility as a function of 
magnetic field at zero temperature was studied in Ref. |43|. 

In this work, we study the general problem of meso- 
scopic fluctuations of thermodynamic properties of the 
grain in the presence of both pairing and exchange cor- 
relations assuming spin-orbit coupling is negligible. The 
quantum Monte Carlo method and Richardson's solution 
mentioned above are computationally intensive, and are 
less practical in calculating the mesoscopic fluctuations 
of thermodynamic properties for which many realizations 
of the grain must be studied. Richardson's solution also 
becomes less tractable at larger values of the pairing gap 
or at higher temperatures, where a very large number of 
energy eigenvalues is required. 

Here we use a more efficient method to calculate the 
heat capacity and spin susceptibility of the grain at finite 
temperature. The exchange interaction is treated exactly 
using a spin-projection method^i 4 ^ and the correspond- 
ing spin-projected partition functions are calculated in 



the SPA+RPA approach. Number-parity effects are cap- 
tured by a number-parity projection. This approach is 
particularly suitable for studying the mesoscopic fluctu- 
ations. 

The outline of the paper is as follows. In Sec. [TT1 we 
discuss the universal Hamiltonian and briefly review the 
conditions of its validity. In Sec. IIII1 we discuss the cal- 
culation of the canonical partition function, and use it 
to evaluate the heat capacity and spin susceptibility of a 
system described by the universal Hamiltonian. We also 
discuss the stability of the RPA, which is unstable below 
a certain critical value of the temperature. In Sec. lIVl we 
present our results and discuss their physical significance. 
We conclude in Sec. E] 

II. MODEL 

In a chaotic grain, the statistical fluctuations of single- 
particle energies and wave functions follow RMT^ at en- 
ergy scales below Etyi ■ In the absence of spin-orbit scat- 
tering and orbital magnetic field, the relevant ensem- 
ble is the Gaussian orthogonal ensemble (GOE). In gen- 
eral, the matrix elements of the electron-electron inter- 
action in the basis of eigenstates of the non-interacting 
part of the Hamiltonian have a complex structure that 
depends on the particular realization of the one-body 
Hamiltonian. In the limit of a large Thouless conduc- 
tance <7Th = Et^/S 3> 1, these matrix elements can be 
decomposed into an average and fluctuating partsi 7 ^ The 
average interaction terms together with the one-body 
Hamiltonian are referred to as the universal Hamilto- 
niani 7 -^ The fluctuating (non-universal) part of the in- 
teraction forms an induced two-body ensemble 4 ^ whose 
matrix elements are suppressed by f/fl'Th- 

For a fixed number of electrons, the universal Hamil- 
tonian has the form 

H= E eiclcia-gPtp-jA*, (I) 

i,cr=f,4. 

where 

p i = T,4t c l, p = Y. c ^' ( 2 ) 

i i 

and S is the total spin. The single-particle levels are 
distributed as the eigenvalues of a GOE random matrix 
with a mean single-particle level spacing S. The uni- 
versal interaction terms are invariant under orthogonal 
transformations of the single-particle basis, allowing us 
to write the one-body part in a diagonal form. If the 
particle number N is not fixed, then the charging term 
E C N 2 should also be included in Eq. ([T]). The bandwidth 
of the model ~ N sp S, where N sp is the number of single- 
particle levels, should satisfy the conditions N sp 5 ^> A 
and N sp S ^> kT at temperature T (where k is the Boltz- 
marni constant). 

The condition of a large Thouless conductance g^h 3> 1 
guarantees that the number of single-particle levels that 
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follow RMT statistics within a Thouless energy window 
around the Fermi energy is sufficiently large, and that 
the nonuniversal correction to the interaction can be ig- 
nored. We assume that the Thouless energy is larger 
than the bandwidth so that RMT is applicable in the 
full model space. This implies the conditions £/rh 3> kT 
and -Erh 3> A. For ballistic grains, the latter assump- 
tion is equivalent to L <C £o , where L is the linear size of 
the grain and £o is the superconducting coherence length. 
This condition also allows us to ignore the spatial fluc- 
tuations of the gap A within the grain, i.e., to treat the 
grain as a zero-dimensional object with respect to the 
fluctuations of the order parameter. 

The pairing constant g in the universal Hamiltonian 
depends on the number of single-particle orbitals in 
the model space. To reduce the computational effort, 
we choose N sp to be smaller than the number of or- 
bitals within the physical window of the Debye frequency 
around the Fermi energy. In doing so we renormalize the 
value of g according to 34 i 46 



stant k = 1) can be written as 
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arcsinh ( 



N ap /2 
A/S 



(3) 



where we have taken the Fermi level to be in the mid- 
dle of the single-particle spectrum (i.e., we assume half- 
filling). Depending on the temperature, we choose JV sp 
to be between 30 and 60 in our calculations to ensure the 
condition N sp S ^> kT. In our studies, A/S < 5.0 so the 
condition N sp S 3> A is also satisfied. 

To study the mesoscopic fluctuations of thermody- 
namic observables, we generate a large number (<~ 1000) 
of realizations of the single-particle spectrum in Eq. (JTJ 
and calculate these observables for each of them. 



III. THEORY 

In this section, we present approximate analytical re- 
sults for the partition function and the spin susceptibility 
of a system described by the Hamiltonian (JTJ for a given 
realization of the single-particle spectrum. First, we show 
that both quantities can be related to the S^-projected 
partition function in the absence of exchange (S z is the 
spin projection along the z axis). Second, we present 
an auxiliary-field path-integral formalism for treating the 
pairing interaction and explain how the integral is eval- 
uated in the SPA+RPA method. Third, we discuss the 
inclusion of number-parity projection. 



Z(J s )=Tre-^ = Tre-^ BCS - J * s 

r 3J,S(S\l) Trri ( : ;//;„ 

S 



^ e /3J s S(S+l) Trs (VMBCs) , (4) 



where Trg is the trace restricted to states with fixed spin 
5 1 , and -Hbcs is the BCS-like pairing Hamiltonian 



HbCS= ^2 e *4r c i<r ~ dP^P- 

t,CT=t,J. 



(5) 



Similarly, the spin susceptibility (at zero external Zeeman 
field) can be written as 



4/3 M | £s S(S + l)e^(s+DTr s (e~^ 
~ 3 Z(J7) 



(6) 



where S z is the spin-component operator along the z di- 
rection and is the Bohr magneton. 

For a scalar operator (i.e., an operator that is rota- 
tional invariant in spin space) 



Tr s X = (25 + 1) (Tr Sz=s X - Tr Si=s+ iX 



(7) 



where Tr s z =m denotes the trace restricted to states with 
a fixed spin component S z — M (while the spin quantum 
number is no longer fixed). Using Eq. ([7]), we can express 
the partition function (j4|) and spin susceptibility ([6]) of 
a system described by the Hamiltonian (JTJ) in terms of 
the Sz-projected partition function of the corresponding 
system in the absence of exchange interaction as 



Z(J s ) = ^2(2S + l)e 



0J B S(S+1) 



(Zs z= s - Zs z =s+i) , 



and 

X = 

Here 



(8) 



3 Z(J S ) 



* + 1)(2S + 1) 



x e 



/3J S S(S+1) 



(Zs z =s ~ Zs z =s+i) ■ (9) 



Zs z 



--M 



Tr 



s z =m e 



-/SHbcs 



(10) 



is the S'z-projected partition function of the BCS-like 
Hamiltonian (|5|). 



A. Spin-projection method 



B. Auxiliary-field path-integral formulation 



The universal Hamiltonian ([T]) commutes with the 
total-spin and particle-number operators. Consequently, 
the corresponding partition function at temperature T = 
1//3 (here and in the following we set the Boltzmann con- 



In the following we consider a partition function of the 
type 



Z\ = Tr>C/ , 



(11) 
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where U is the propagator of the BCS Hamiltonian 

fj = e -p(H BC s-HN) > ( 12 ) 

and Tt\ denotes a trace with certain restrictions (e.g., 
fixed S z ). Using the auxiliary-field path- integral repre- 
sentation of the propagator U, we discuss the static-path 
approximation (SPA) as well as the quantal RPA-like cor- 
rections around each static fluctuation. We have included 
a term in the propagator (Ti~2l) because further calcu- 
lations of the partition function ([TU]) will be carried out 
in the grand-canonical formalism. 

The interaction term in the BCS Hamiltonian ([5]) can 
be decoupled by means of a Hubbard-Stratonovich trans- 
formation! 12 ' 13 In this transformation, we express the 
propagator (|12[) as a functional integral over a complex 
auxiliary field A(r) that depends on imaginary time r 
(0 < r < (3) as follows* 2 ^^ 



~ ~ - f dT\A(r)\ 2 /g „ 

U= I V[A,A*]e o U A 



Here, 



(13) 



The time-dependent part of Eq. (|15l) can be written as 

A*_ rCii CiA . (19) 



V(t) = - J>- 



Since we are interested in the correction to the SPA prop- 
agator (|18p . it is natural to use an interaction represen- 
tation of V{t) with respect to the unperturbed static 
Hamiltonian H/\ 

tW^e^^e-^o . (20) 

To obtain the expression for Vl n t, it is convenient to 
work in the Ao-dependent basis that diagonalizes the 
static Hamiltonian H a . We will refer to this basis as the 
quasiparticle representation (even though this is strictly 
the case only for the saddle-point value of Ao). For a 
given value of Ao, the quasiparticle operators are re- 
lated to the original particle operators via the Bogoliubov 
transformation 



c , it = ma^ + Vi e 10 , 
c-4 = UiCtii - Vi e ie a| f . 



(21) 



I e 



(14) 



is the propagator for the one-body Hamiltonian 



n 



' Ac lt c 4- A*c l4 c iT + | 



(15) 



and T denotes time ordering. The measure 2? [A, A*] is 
defined to preserve the normalization (i.e., U = 1 when 

D[A, A* 



M „ 

lim TT / 

m=l J 



2ng 



dA m dA* m , 



(16) 



where Af3 — jj. 



The auxiliary field A(r) can be separated into its static 
and r-dependent parts by a Fourier series 



A(r) = A + ^A r e^, 

r=iO 



(17) 



where ui r = 2nr/ '(3 (r integer) are the bosonic Matsubara 
frequencies. In the SPA, A(t) is replaced by Ao in the 
Hamiltonian (p~5]) . and the propagator (fl4l is approxi- 
mated by 



Ua 



(18) 



where Ha is the Hamiltonian (TT51 for a static field A = 
A . The BCS theory can be derived by applying the 
saddle-point approximation to the SPA integral. 



Here 



-(1-7,), = argA o , (22) 



where 



7i 



u - £ 

fj, 2 



E, 



and Ei are the quasiparticle energies 

Ej = \ / | ( ; — fX ■ 



Ac I 



(23) 



(24) 



In this quasiparticle basis, Ha an d ^( T ) have the forms 



#A 

and 
V(r) 
. |Ao| 



y^lej - fj, - Ej + Ej (a^a^ + a^a.j 



(25) 



2 ^ 



Ei 
with 



A ir a| t at ; +A*_ r a 2i a lt 



1 - 4^ - at f a it 



A 4r = (7i + l)A r + (ji - l)e 2i6 A 



(26) 



(27) 



The dependence of V(t) on 9 can be eliminated by an 
appropriate gauge transformation of A r and the quasi- 
particle operators. Therefore, the integrand in the SPA 
integral depends only on the absolute value of Ao. 
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The form of Vi n t(r) in the quasiparticle basis is the 
same as of V(t) in Eq. (|26[) with a and replaced by 

a la (r) = a ta e~ E ^ , <4(r) = a\ a e E ^ . (28) 
In the interaction representation, 

U A = U Ao U int , (29) 

where 



-fdrV lnt (r) 

U lnt = Te o 



(30) 



As a result, the partition function can be written as 

oo 

Z x = ^ [ d |A |Vf |Ao ' 2 Z X (A ) C X (A ) , (31) 
o 



where 



Z A (Ao) - Tr x U Ao 



(32) 



is the partition function for a static value Ao of the pair- 
ing field, and 

is the quantum correction. The average in Eq. ([33]) is 
defined with respect to the static-field propagator 



A 



Tr. 



(y Ao a) 



a,a Tr x lI Ao 
and the integration measure is 

/3dA r dA* r 



V'[A r ,A* r ] = l[ 



2tt 9 



(34) 



(35) 



Equation (1311) for the partition function Z x is ex- 
act. In the SPA, C\ = 1 or, equivalently, V"(t) = 0. 
The SPA+RPA approximation is obtained by evaluating 
the integral over A r and A* in Eq. (|33[) in the saddle- 
point approximation assuming small-amplitude fluctua- 
tions. To this end, we write 

ln/^nt) « I f f dTdT'(TV int (T)V int (T')) X . Ao , 

\ / A,Ao 2 J Q J Q 

(36) 

where we have assumed that J dT(Vi n t(r))x, Ao = 0. 
This is valid provided the projection in Tr x conserves the 
quasiparticle occupation number, which is the case for 
the Sz -projected trace but not for the canonical (particle- 
projected) trace. The RPA correction factor (|3"3"f is then 
given by 

C* PA (A ) = / V'[A r , A*] e -4£^o 



xexp U J P J^Tdr' (TV u ,ir)\l,',r' ) 



A,A 



(37) 



Higher-order corrections can be obtained by including 
more terms in the cumulant expansion (|36j) . 



C. Canonical and number-parity projections 

The canonical partition function can be related to 
the grand-canonical partition by particle-number projec- 
tion^ However, an exact particle-number projection "in- 
side" the integral (|31|) for each value of the static field 
Ao leads to a complicated expression since the particle- 
number operator N does not commute with the static 
Hamiltonian H Aq £Z. 

In the following we will carry out the particle-number 
projection in the saddle-point approximation. However, 
to account for important odd-even effects, we project on 
the number parity of electrons. The trace in (TTT|) will be 
restricted to S z and the number-parity quantum number 
77(77 = +1 for an even number of particles, and 77 = — 1 
for an odd number), i.e., A = 77, S z . 

The canonical partition function for TV particles can be 
obtained from Eq. (fTTj) by particle-number projection 



Zn.x = 



da 
27 



-iaN 



fidfj, j fldJApj e -0lF x (u,Ao)+nN] Cx (Ao) , 
2ni J g 




where 



e -/J^(M.Ao) = e -(/3/ 9 )|A | 2 Za(Ao ). 



(38) 



(39) 



Here Z X (A ) and C X (A ) are calculated for a chemical 
potential ix. We denote by F the grand-canonical free 
energy (without the restriction A) given by 



-PF(n,Ao) = e ~(l3/g)\Ao\ 2 TTe -l3H Ao 



fJ^-^-^cosh 2 ^. (40) 



= e -(/3/ S )iA r 



We exchange the order of integrations in Eq. (I38|) and 
for each value of A evaluate the integral over /i by the 
saddle-point approximation. The latter is applied to the 
function e~^^ F ^ :A( ^ +fJ/N ^ considering the remaining part, 
e~P( Fx ~ F )C\(Aa), as a pre-factor that does not affect the 
saddle-point integration. We then obtain^ 



oc 


0d|A 


-/ 




9 


/2vr 


8 2 F s 


\J 


d[i 2 / 



-1/2 



^[f a ( m ,Ao)+mA1 Ca ( Ao ) ; ( 41 ) 
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where 



c^F = _ ^ /3E t { e t - n - § ) 2 + |A | 2 smh(/3^) 



2j B?co S h 2 (^) 



(42) 

and ji is a Ao-dependent chemical potential determined 
by 



The number-parity projection is carried out using the 
projector 



^ = o 1 + V? 



tN 



(44) 



where r\ = ±1, depending on the number-parity of elec- 
trons. If this operator is inserted in the grand-canonical 
trace, only states with even (odd) number of particles 
will be taken into account for rj = 1 (77 = — 1). 



D. Number-parity and S^-projected static partition 
function 



Here we discuss the calculation of the static partition 
function (|32j) . when A corresponds to the projections on 
number parity and S z , i.e., A = 77, S z . 

The trace over states with fixed S z can be calculated 
exactly through a discrete Fourier transform as long as 
the maximal value of S z is finite. The corresponding 
projection operator is 



Ps 



1 



25„ 



in— — o ma > 



(45) 



where SWiax is the maximal possible value of the spin and 
<j) m = 2"Km/(2S max + 1) are quadrature points. We use 
this operator in accordance with the number-parity pro- 
jection, i.e., the values of m are integers or half-integers 
for even or odd number of electrons, respectively. Be- 
cause our goal is to obtain an expression for the canon- 
ical partition function, the value of S max is determined 
by the number of particles, rather than by the size of the 
Hilbert space. 

The spin projection operator S z commutes with H& , 
so we can write 



Tr 



(V 



Tr(P„e^ ffA ° + ^ m5 * 



(46) 



The second-quantized forms of S z and P v remain in- 
variant under Bogoliubov transformation (|2"Tj) . Conse- 
quently, the projected partition function for a static pair- 
ing field A is given by 



Z„ A (A ) = £ 



25*77 



(47) 



where 

^(°.*m)(A ) = Tr ( e -^n+^^ 



1 + e 



, (48) 



and 



(A ) = Tr (. 



e iirN e -f3H &Q +i4> m S s 



-P(ei-fi-Ei) 



(49) 



The last two expressions were obtained using the grand- 
canonical formalism. 



E. The RPA correction 

The RPA correction factor ([57| is calculated in Ap- 
pendix |XJ For A = 77, S z , it is given by 



Cf A (Ao) = n[det^K)r x 



(50) 



r>0 




= l dlnZ A (A ) 



(51) 
(52) 



Here w r = 2irr/f3 are the positive bosonic Matsubara 
frequencies, and 7$ are given by Eq. (|23l) . Each factor in 
the product ([50]) is obtained after a Gaussian integration 
over A r , A*, A_ r and Al r . This integral converges if and 
only if the real parts of both eigenvalues of the matrix 
A{ui r ) are positive. Since A{ui r ) is a 2 x 2 real matrix, 
this is equivalent to 



det^4(w r ) > and tvA(uj r ) > . 



(53) 



For a given Ao, the RPA breaks down below a certain 
temperature for which one of the conditions in Eq. (|53[) is 
not satisfied. It is then necessary to include higher-order 
cumulants in the expansion (1361) . The highest tempera- 
ture for which the RPA breaks down for at least one value 
of Ao is known as the SPA+RPA critical temperature T* . 
The SPA+RPA approach is thus valid for temperatures 
above T». Numerical simulations show that T* increases 
with the BCS gap A and becomes of the order of S for 
A/ 5 ~ 5. It has been recently proposed^ that the above 
instability can be avoided by treating non-perturbatively 
a low-energy collective mode. 



In Appendix [B] we show that the RPA correction fac- 
tor (|5H)) can be written as 

where ±f2j are the eigenvalues of the 27V sp x 2N sp RPA 
matrix (-/V sp is the number of single-particle orbitals) 

/ 2EA - ZM-Yili + 1) -|/Ai(7i7j - 1) A 
^ §/A<(7i7i - 1) §/Ai(7i7j + 1) - y 1 ■ 

(55) 

Equation (|54[) is more practical since it does not contain 
an infinite product and we have used it in our calculations 
below. 

Equations ([50]) and (|54l) are valid not only for the 
restricted partition function with A = T],S Z , but 
also for the grand-canonical partition functio n 25 i 26 ' 30 and 
for the number-parity-projected partition function^ In 
all of these cases, the correct Z\(A ) must be used in 
Eq. ([52")) to define fxi- These expressions, however, are 
not applicable for the canonical projection. 

F. Summary 

In summary, we use Eqs. © and © to express the 
iV-particle partition function Zn and spin susceptibility 
X of a system described by the universal Hamiltonian ([T]) 
in terms of the number-parity and ^-projected partition 
function Z N ^ tSz of a system described by the BCS-like 
pairing Hamiltonian ([5]) [note that in Eqs. (JS) and (j9|) we 
replace Z$ z by Zat^sJ- The partition function -Zjv,?7,Sj 
is calculated in the SPA+RPA using 

7 /3d|A | 2 {2ir d 2 F\- 1/2 

Zw ° ~ J [J w ) 

o 

x e -(/V 9 )|A | 2 e -P»»z v ,s.(*o) ^(Ao) • (56) 

Here Ao denotes a static fluctuation of the order param- 
eter, and r\ is the number parity (77 = 1 for even N and 
— 1 for odd N). The partition function Z v> s z (^o) and 
the RPA correction C3^(A ) are given by Eqs. (03) and 
(|50"| or (j54"|) . respectively. The second partial derivative 
d 2 F/d^i 2 is given by Eq. (J32), and the chemical potential 
fi for a given static fluctuation Ao is determined from 
Eq. (|43l) . The heat capacity C is obtained numerically 
from the partition function Zn (T) as a function of tem- 
perature. 

IV. RESULTS AND DISCUSSION 
A. Accuracy of the method 

We first discuss the accuracy of the number-parity pro- 
jected SPA and SPA+RPA methods. To this end, we 
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A/8 = 3.0, J /8 = 0.5 
15i 1 1 1 1 1 1 1 1 1 




FIG. 1. Heat capacity C vs T/S for even and odd grains 
with A/S — 3.0, J s /S — 0.5, and for a specific realization 
of the GOE single-particle spectrum in ([T|l. The results cal- 
culated in the number-parity projected SPA (open symbols) 
and SPA+RPA (solid symbols) for even (circles) and odd 
(squares) grains are compared with exact canonical results ob- 
tained by Richardson's solution (solid line for the even grain 
and dashed line for the odd grain). The results based on 
Richardson's solution use all eigenvalues below a cutoff of 
~ 30 8 and are no longer accurate for temperatures above 
T ~ 1.55. 



have used the exact solution of the Hamiltonian ([TJ , mod- 
ifying Richardson's solution for the BCS-like Hamilto- 
nia n 35 i 36 to include the exchange interactionJS The num- 
ber of many-body eigenstates that contribute to the par- 
tition function increases rapidly with temperature, and 
so does the required computational effort. In practice, 
we compute only the energy eigenvalues below a cutoff 
of ~ 30 6. These exact calculations are then accurate at 
sufficiently low temperatures where the contribution of 
eigenstates with energy above 30 5 is negligible. 

The comparison between the exact and approximate 
calculations is demonstrated for a given realization of the 
single-particle spectrum in Figs. [TJ and [2] for the heat 
capacity and spin susceptibility, respectively. We show 
results for both even and odd number of electrons. 

We observe that the number-parity projected 
SPA+RPA (solid symbols) improves significantly the 
number-parity projected SPA (open symbols) and 
provides accurate results for both the heat capacity and 
spin susceptibility. In particular, the RPA correction 
is important for the spin susceptibility at larger values 
(closer to 1) of the exchange coupling J s /S. The example 
shown on Fig. [2] demonstrates that the number-parity 
projected SPA results can be qualitatively wrong, 
whereas the inclusion of the RPA correction factor gives 
much more accurate results. 
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A/8 = 0.5, J /8 = 0.5 



A/5 = 0.5 



A/5 = 1.0 



A/5 = 3.0 




FIG. 2. As in Fig. [TJbut for the spin susceptibility \ normal- 
ized by the Pauli susceptibility \p — ^A/S vs T/8 and for a 
grain with A/5 = 0.5, J 3 /S — 0.5. Symbols and lines follow 
the same convention as in Fig. [1] 




B. Heat capacity 

There are two major number-parity-dependent signa- 
tures of pairing correlations in the heat capacity: the 
heat capacity for an even particle number is suppressed 
at low temperatures and enhanced at intermediate tem- 
peratures when compared with the heat capacity for an 
odd particle number (see Richardon's solution results in 
Fig. [JJ) . The low-temperature effect is not accessible by 
the method we are using because the RPA becomes un- 
stable at low temperatures, and we concentrate below on 
number-parity effects in the intermediate temperature re- 
gion, ft is known that, in the absence of exchange, the 
characteristic temperature of this region is determined 
by the scale that is the largest between S and Ai2i In the 
BCS regime (i.e., large A/S), this effect occurs around 
the BCS critical temperature, while in the fluctuation- 
dominated regime A/S < 1 it occurs at temperatures 
higher than the BCS critical temperature. The even-odd 
effect becomes more prominent when the size of the grain 
and consequently A/ 5 increase. The heat capacity for 
even particle number starts displaying a shoulder around 
A/S ss 3.0, which eventually develops into a sharp peak 
in the bulk limit A/S 3> 1. Here we investigate how this 
picture is affected by a non-zero exchange interaction and 
by mesoscopic fluctuations. 

The results for the heat capacity are shown versus T/S 
and both number parities in Fig. [3]for A/S = 0.5, 1.0, 3.0 
and J s /S = 0,0.4,0.8. The symbols and vertical bars 
are average values C and standard deviations SC, re- 
spectively, calculated from an ensemble of 1000 random 
matrices describing the one-body part of the Hamilto- 
nian (TTJ). The lines (solid for even and dashed for odd 
number of electrons) are obtained for the equally spaced 
single-particle spectrum in the Hamiltonian ([1} with level 



FIG. 3. The heat capacity C vs temperature T/5 for an even 
grain (solid lines, circles) and for an odd grain (dashed lines, 
squares). Results are shown for grains with A/S — 0.5 (left 
column), A/8 = 1.0 (middle column), and A/6 = 3.0 (right 
column) and with J s /8 — (top row), J s /S = 0.4 (middle 
row), and J s /8 = 0.8 (bottom row). The symbols and verti- 
cal bars describe, respectively, the average value C and stan- 
dard deviation SC of the heat capacity (where an ensemble of 
single-particle spectra are sampled from the GOE) . The lines 
correspond to an equally spaced single-particle spectrum in 
the Hamiltonian (TTJ and the dash-dotted lines are the grand- 
canonical BCS results (where applicable). 



spacing 5. Figure 2] shows the standard deviation SC ver- 
sus T/S for the same cases as in Fig. [3l 

We observe that the exchange interaction can suppress 
the odd-even effects in the heat capacity and shift them 
to lower temperatures. This is particularly evident if 
A/S < 1. For A/S = 3.0, a higher value of the exchange 
coupling constant is required to make a visible change. 
Even for J s /S = 0.8 (which is close to the Stoner insta- 
bility threshold), the number-parity effect is shifted to 
lower temperatures slightly. Only the right side of the 
even number-parity shoulder is suppressed, whereas the 
left side is not. As a result, the shoulder transforms into 
a peak. 

This behavior is consistent with the ground-state phase 
diagram 40 of the grain, according to which the ground 
state for an even particle number is fully paired for small 
J s /S and the value of J s /S required to polarize the grain 
increases with A/S. For A/S = 3.0, this value of J s /S 
is close to the Stoner instability threshold. For smaller 
values of J s /S, the excited states with non-zero spin be- 
come important only at sufficiently high temperatures. 
At lower temperatures, the dominant contribution to the 
heat capacity comes from the zero spin levels whose en- 
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T/8 

FIG. 4. The standard deviation SC of the heat capacity 
(shown by vertical bars in Fig. [3]) versus temperature T/8 
for an even grain (solid lines) and for an odd grain (dashed 
lines) with the same values of A/5 and J s /S as in Fig. [3] 



ergy is independent of J s /5. This is consistent with the 
weak dependence on J s /8 of the left side of the even-case 
shoulder for A/8 = 3.0. At temperatures that corre- 
spond to the right side of the even-case shoulder, non- 
zero spin configurations are more important and lead to 
visible change in the heat capacity. For A/ 5 < 1, the 
excitation energies of the states with non-zero spin are 
lower and the heat capacity is more sensitive to exchange 
correlations at lower temperatures. 

In an ultra-small grain, the mesoscopic fluctuations 
of observables can be significant. For example, if an 
odd-even signature of pairing correlations is studied by 
carrying out measurements in many samples with un- 
known number parity and then determining the distribu- 
tion of the observable, such a number-parity effect may 
be washed out when the fluctuations are large. As can 
be seen from our results, this can happen if A/ 6 is suf- 
ficiently small or J s /6 is sufficiently large. In the ab- 
sence of exchange, this occurs at A/ 8 < 0.5, and, for 
J s /8 ~ 0.5, number- parity effects already disappear be- 
low A/ 6 ~ 1. However, for A/ 5 = 3.0, the odd-even 
effect is not suppressed by the fluctuations even at rela- 
tively large J s /8, 

As can be seen in Fig. |4] the mesoscopic fluctuations of 
the heat capacity at high temperatures are only weakly 
dependent on the pairing gap, exchange coupling con- 
stant and the number parity. At intermediate temper- 
atures, when the heat capacity is enhanced for an even 
particle number, the fluctuations in the even case are 
stronger than in the odd case and are characterized by 
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FIG. 5. The spin susceptibility % m the units of the Pauli 
susceptibility \p = 2/is/<5 vs temperature T/6 for even and 
odd grains with A/5 = 0.5 (left column), A/8 — 1 (middle 
column) and A/5 = 3.0 (right column) and with J s /S — 
(top row), Js/S = 0.3 (middle row) and J s /S = 0.6 (bottom 
row). Symbols and lines follow the same convention as in 
Fig. [3] but for the spin susceptibility. 

a peak. The position and height of this peak are almost 
independent oi A/8 for A/ 8 < 1, but increase with A/ 8 
for A/8 > 1. In the presence of exchange correlations, 
the peak is shifted to lower temperatures. This is consis- 
tent with a similar shift of the odd-even signature in the 
heat capacity. In addition the size of the peak increases 
with J s for A/ 8 > 1. 



C. Spin susceptibility 

It is known that pairing correlations (in the absence of 
exchange) suppress the spin susceptibility for both num- 
ber parities. For an odd particle number, this suppression 
together with the low-temperature Curie-like divergence 
(~ leads to a re-entrant behavior. For an even 

particle number, the spin susceptibility increases mono- 
tonically with temperature. Exchange correlations and 
mesoscopic fluctuations may affect this behavior. 

We express the spin susceptibility \ m the units of the 
Pauli susceptibility xp = 2/^/5, where \xb is the Bohr 
magneton. For given values of A/8 and J s /8, the ratio 
x/xp is expected to be a universal function of T/8. In 
the high-temperature limit, x/xp does not depend on 
A/ 8. For an equally-spaced single-particle spectrum and 
in the absence of exchange, it approaches 1 in that limit. 

The results for the spin susceptibility versus T/8 are 
shown (for both number parities) in Fig. [5] for A/ 8 = 



in 
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FIG. 6. The standard deviation 8x1 Xp of the spin suscepti- 
bility (shown by vertical bars in Fig. [5| vs temperature T/8 
for an even grain (solid lines) and for an odd grain (dashed 
lines) with the same values of A/ 8 and J a /8 as in Fig. [5] 



0.5, 1.0, 3.0, and J s /8 = 0, 0.3, 0.6. Symbols and lines 
follow the same convention as in Fig. [3] The standard 
deviation Sx/xp is shown in Fig. [6] versus T/8. 

The most visible effect of the exchange interaction is 
the enhancement of x/xp as higher spin states are shifted 
down in energy. Exchange correlations can also shift the 
odd-even effects to lower temperatures and even elimi- 
nate the odd-case re-entrant behavior for A/<5 < 1. How- 
ever, at larger values of A/ 8, exchange enhances the re- 
entrant behavior. This effect is similar to what we ob- 
served for the even-case heat capacity where a shoulder 
changes into a peak for larger values of A/ 8 (see the right 
column of Fig. [3|). 

The mesoscopic fluctuations of x/xp increase with 
decreasing A/5 or increasing J s /5. In the fluctuation- 
dominated regime A/ 5 < 1, they can become especially 
strong at larger values of J s /5 and can hinder the obser- 
vation of odd-even effects. When compared to the heat 
capacity results, higher values of J s /5 or smaller values 
of A/ 5 are required to hinder the odd-even effects. 

The large mesoscopic fluctuations of the spin suscep- 
tibility for A/ 5 < 1 and large values of J s /5 may be 
attributed to the large dispersion of the magnetization 
of the system. This is confirmed by studying spectra 
of individual samples in that regime using Richardon's 
method. Examples are shown in Fig. [7J Large spin sus- 
ceptibility values are obtained in samples in which the 
excitation energies of states whose spin is different from 
the ground-state spin are particularly low. The proba- 
bility to have such a sample is enhanced when pairing 
correlations are weaker and/or exchange correlations are 
larger. 

At relatively large values of A/ 5, the fluctuations in- 
crease monotonically with temperature (see right column 
of Fig. [6]). However, for A/ 5 < 1, the fluctuations in the 
spin susceptibility have a maximum at smaller values of 
Js/5. It is not clear if this is the case at large values J s /5 
because we cannot access very low temperatures by our 
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FIG. 7. The spin susceptibility x/Xp f° r an even grain with 
A/8 = 0.5 and J s /8 — 0.5 for three different RMT realizations 
of the single-particle spectrum in The table shows (for 
each of the three realizations) the lowest excitation energy 
Es / 8 for a given spin S. 



method. However, we expect the fluctuations to remain 
strong in the limit T — > in the odd case because of 
the Curie-like behavior ~ Sq/T, where So is fluctuating 
ground-state spin. 

We note that the equally spaced single-particle spec- 
trum does not always describe the average behavior of the 
spin susceptibility. In the regime where the mesoscopic 
fluctuations are strong, the spin susceptibility calculated 
for the equally-spaced spectrum (solid lines in Fig. [5]) is 
smaller than the spin susceptibility obtained by averag- 
ing over all samples (circles) and may be qualitatively 
different (see, e.g., the case A/8 = 0.5 and J s /8 = 0.6 in 
Fig. EJ. 



V. CONCLUSIONS 

In conclusion, we have studied the thermodynamic 
properties of ultra-small chaotic metallic grains with a 
large dimensionless Thouless conductance in the presence 
of both superconducting and ferromagnetic correlations. 
We have used the so-called universal Hamiltonian (JXJ) as 
our model, in which the one-body part is sample-specific 
and modeled by RMT, while the dominating interaction 
terms are universal. Sample-to-sample fluctuations of 
the interaction are suppressed and ignored in the limit of 
large dimensionless Thouless conductance. The exchange 
interaction has been treated exactly by means of a spin- 
projection method, while the pairing interaction has been 
treated in a path-integral approach in which all static 
fluctuations of the pairing gap and small-amplitude time- 
dependent fluctuations around each static value of the 
gap are included (SPA+RPA method). Particle- number 
projection is approximated in the saddle-point approxi- 
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mation while number-parity effects are preserved using 
an exact number-parity projection. The method is effi- 
cient and very accurate when compared to exact canoni- 
cal calculations. However, it cannot be used at very low 
temperatures, when the RPA correction becomes unsta- 
ble. This limitation can potentially be overcome using 
the method developed in Ref. H^. 

We have found that the exchange interaction shifts the 
number-parity-dependent signatures of pairing correla- 
tions (such as the enhancement of heat capacity in the 
even grain and the re-entrant behavior of the spin sus- 
ceptibility in the odd grain) to lower temperatures. In 
the fluctuation-dominated regime A/ 5 < 1, these signa- 
tures are suppressed by exchange correlations. However, 
at sufficiently large values of A/ 5, exchange correlations 
have the opposite effect, i.e., the heat capacity of the even 
grain develops a peak, and the re-entrant behavior of the 
spin susceptibility in the odd grain is enhanced. 

Mesoscopic fluctuations of thermodynamic observables 
can further hinder the odd-even effects for sufficiently 
small A/ 8 and large J s /S. The mesoscopic fluctua- 
tions of the spin susceptibility are especially large in the 
fluctuation-dominated regime A/S < 1 for values of J s /S 
above ~ 0.5 because of the large dispersion of the mag- 
netization. 

It would be interesting to extend our work to the study 
of granular metals^ i.e., arrays of metallic nanoparticles 
that are coupled via tunnel junctions. For weakly cou- 
pled grains and when the charging energy Ec satisfies 
Ec A, the majority of grains are in the Coulomb- 
blockade regime kT <C Ec with suppressed inter-grain 
tunneling. These Coulomb-blockaded grains provide the 
dominant contribution to the thermodynamic properties 
of the granular metal at low temperatures. Therefore, the 
values of thermodynamic observables of a granular metal 
(per grain) at kT <C Ec can be effectively calculated by 
averaging the observables of individual grains over differ- 
ent random-matrix realizations and the number parity of 
electrons. We note that number-parity effects must still 
be taken into account since they lead to effects that could 
be missed in grand-canonical calculations. An example 
is the Curie-like divergence in the average spin suscepti- 
bility. 
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Appendix A: Calculation of the RPA correction 
factor 

In this appendix we calculate the RPA correction fac- 
tor (|37|) with A denoting simultaneous projections on the 
number parity r\ and spin component S z . The corre- 
sponding projection operators, given by Eqs. (|44[) and 
(|4"S"|) . respectively, have the same form in both single- 
particle and A -dependent quasiparticle representations. 
This indicates that they commute with the quasiparti- 
cle occupation number operator a ia ai a for each i and a. 
Therefore, if Vint(i~) is written in the quasiparticle repre- 
sentation according to Eq. (|26|) with r-dependent oper- 
ators, non-zero contributions to the correlator of V[ n t (t) 
in Eq. (pTT)) are possible only from the product of ctj-t a i| 
and a^aj-j- and from the product of two terms of the form 

1 - <4j.°U ~ a lf a «t taken from both Mnt(r) and V^t'). 
In the latter case, the two terms are r-independent and 
identical. Thus time-ordering can be omitted, and inte- 
gration over t vanishes because of the e luirT factor. Con- 
sequently, we obtain 

@ & 

\ j J 'drdr' (TtWT)lW/)) a 
o o 

(Al) 

where 

B i (Ty) = (Ta^(T) a l(T)a il (T')a it (r'))^^ . (A2) 

Wick's theorem cannot be applied directly to the cor- 
relator Bi(r, t') because of the projections. To proceed, 
we use the following two identities 

( a 4^t a I t a !;) A Ao = e2/3jEi ( a lAi a *i a *t) x Ao ( A3 ) 

and 

(l- e 2 ^)(«l t <^«it) AAo 

= -l + («W) A)Ao + («W) AAo - (A4) 

The second identity follows directly from the first one 
and the anti-commutation relations. To derive the first 
identity, we write the projected average (j34|) of an ob- 
servable A at a given static field in the form 

where C$ x are certain coefficients and 

U) = — ^ pf. (A6) 
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The sum in Eq. (|A5I) is over quadrature points <f)\ suit- 
able for the projection A. In deriving (|A5|) . we have 
used expressions (|44|) and (|45|) for the projection op- 
erators and replaced e mN by e 2mS * (N and 2S 2 have 
the same parity). The expectation value in Eq. (|A6j) 
is grand canonical with respect to the one-body Hamil- 
tonian Ha — i4>\S z /l3 (Ha commutes with S z ), and 
therefore can be calculated using Wick's theorem. We 
find 



04^0^) = (1 - nf t A )(l - n$) 



e 2 ^(4AlW)^ ■ ( A7 ) 



where 



a\„a ic 



A + 1 



(A8) 



Using Eqs. (|A5|) and (|A7|) . we obtain the relation (|A3|) . 

We can now evaluate the integrals on the right-hand 
side of Eq. dXTJ with the help of Eq. d2SJ) and the iden- 
tities (El) and (HI to find 



& 







r' T 'Bi(T,T') = ^al^a^ 



AA 















a it a a a ii a it 



AA 



dre (i Wr +2J3 i )r 



(iui r ,+2Ei)T _ j 2 _ e 2f)E le -(iu r .,+2E i )T 



iu> r i + 2Ei 



icu r r + 2Ei 



p6 r 



iu> r + 2E. 



■hi, (A9) 



where 



^ = i-((«W+«W)) AAo = ^ 



191nZ A (A ) 



dEi 



(A10) 

parts of A(r)e~ l& , respectively, we change variables 

A r = e t6 (<7 lr + ia 2r ) , (All) 
where <Jkr is the Fourier transform of crfe(r) (<7fc,-r = 



<xL.). The quantity A ir in Eq. (1271) can then be written 
as 



A,„. = 2e i8 (jtair + ia 2r ) 



(A12) 



Denoting by oi(r) and ct 2 (t) the real and imaginary an d the integration measure as 



r>0 



n 9 



(A13) 

Consequently, the RPA correction factor is given by 



C^ A (A ) = / V'[a lr ,a 2r ]exp 



"(2/3/5) £ 



r>0 



9 



.9 



4- 4K 2 



0"2r 



4E 2 ^^ ( a ^(T 2r - a lr a 2r ) 



= J[ [dct A(a 



(A14) 



r>0 



where the matrix A(w r ) is defined in (f5Tj) . Note that 
in general Eq. ([50|) is valid as long as Eq. (|A3|) holds or 
the projection operator in Z\ conserves the quasiparticle 
occupation numbers. 



Appendix B: Relation of the RPA correction factor 
to the RPA matrix 



As a function of uj, det^4(w) in (|A14[) has poles at 
uj = ±2iEi, which could be either first or second order. 
A second-order pole at ±2iEi can only arise from prod- 
ucts of two terms that contribute to the matrix elements 
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of A(u>) and have (AE 2 +u> 2 ) in their denominators. How- 
ever, when the sum of these products is written as a ratio 
of polynomials, a partial cancellation between denomina- 
tor and numerator results in first-order poles. 

Since all the poles of det A(ui) are first order, it can be 
written as 



det A(u) 



AE 2 ) 



2\ ' 



(Bl) 



where P(u>) is a polynomial of degree 2iV sp (N sp us the 
number of single-particle orbitals). The roots of P(ui) 
and det A(oj) come in pairs ±if2j (det A(uS) is a function 
of w 2 ), and the leading coefficient of P(ui) is equal to one. 
Thus P(u) = l\ t (uj 2 + Q 2 ), and 



detA( Wr )=TT-^ + ^ 



(B2) 



9 \ - 

i 

9 ST 



sinha; = ^Ilr>o(l + x 2 /n 2 r 2 ), we obtain Eq. ([S3| for 
the RPA correction factor. 

Next, we show that ±0; are the eigenvalues of the 
2N sp x 2jV sp RPA matrix ([55]) . Indeed, considering one 
of the eigenvalues fi of this matrix and denoting its cor- 
responding eigenvector by (xu X2i) T , we have 

9 hi 



Xu 



X2i 



2 2Ei - Q 



9 f. 



2 2Ei + Q 



Defining 



3 

E ((7i7j - + (7<7j + l )Xv) ■ 

(B3) 



. u>* + 4Ef 

Using (IA14[) and the infinite product representation we obtain 
I 



»?+ = E^' +*2j)7j 

= E^ ~ X2j) , 



hiV+ +n-)+ o F : (7^7+ -»?-)J = </ 2^ 



2£; - 
/a 



2^ + n 



2E, - ft 



(7i»?++»7-) 



2E, + n 



(7*»7+ - 



'E 



4s 2 - n 2 
ae 2 - n- 



V+ + 9 E 



■ V+ + 9 



E 



AE 2 - ft 2 ?? " 



2Ejfxi 

ae 2 - n 2 



(B4) 



(B5) 



Equation (|B5[) can be rewritten in the form 

B(Q)(jj+)=0, (B6) 

where is a 2 x 2 matrix satisfying det-B(fi) = 

det A(ifl). In general, (^ + ) ^ and Eq. (jB6]) im- 
plies detB(fl) = 0. Therefore, w — ifl is a root of 
det Since the eigenvalues of the RPA matrix (|55|) 

come in pairs ±f2 and their number is equal to the num- 
ber of roots of det A(uj), we conclude that all of the 
in Eq. (|54|) are eigenvalues of the matrix (1551 and vice 
versa. 

The approximation used to calculate the correction 
factor (|5U|) breaks down when at least one of the ma- 
trices A(ui r ) has an eigenvalue with negative real part 
and the corresponding Gaussian integral diverges. As 
discussed in Sec. IIIIE1 the necessary and sufficient con- 
ditions for the SPA+RPA to be applicable are that all 
matrices A(u) r ) have positive traces and determinants 



[see Eqs. (l53|) ]. It is clear from (|B2[) that det A(uj r ) can 
be made negative only when at least one of the fli is 
complex. Moreover, such Qi must be purely imaginary; 
otherwise, ft* also appears in the product of Eq. (|B2[) 
making the combined contribution from 17^ and f2* pos- 
itive. When Qi is purely imaginary, its complex conju- 
gate belongs to the same "pair" ±f2j of eigenvalues of the 
RPA matrix and does not give a separate contribution to 
the product in Eq. (|B2p . As the temperature decreases 
and an RPA frequency becomes purely imaginary, the 
first determinant that becomes singular is the one with 
the smallest u> r , i.e., wi = 2ttT. Consequently, all de- 
terminants det A(uj r ) are positive if |fij| < 2ixT for all 
purely imaginary RPA frequencies f2,; . There is a critical 
temperature T* below which this condition is no longer 
satisfied for some value of the static gap Ao- 

In principle, the condition |f2j| < 2irT is necessary but 
not sufficient. Simulations show that for a small static 
field Ao, the eigenvalues of A(ui r ) may form a complex- 
conjugate pair such that detA(u; r ) positive, while their 
real parts may be negative. This can happen at temper- 
atures of the order of or lower than T». 



1 D. C. Ralph, C. T. Black, and M. Tinkham, Phys. Rev. 2 C. T. Black, D. C. Ralph, and M. Tinkham, Phys. Rev. 
Lett. 74, 3241 (1995). 



14 



Lett. 76, 688 (1996). 

3 D. C. Ralph, C. T. Black, and M. Tinkham, Phys. Rev. 
Lett. 78, 4087 (1997). 

4 F. Kuemmeth, K. I. Bolotin, S. F. Shi, and D. C. Ralph, 
Nano Lett. 8, 4506 (2008). 

5 Y. Alhassid, Rev. Mod. Phys. 72, 895 (2000). 

6 M. L. Mehta, Random Matrices, 2nd ed. (Academic, New 
York, 1991). 

7 I. L. Kurland, I. L. Aleiner, and B. L. Altshuler, Phys. 
Rev. B 62, 14886 (2000). 

8 I. L. Aleiner, P. W. Brouwer, and L. I. Glazman, Phys. 
Rep. 358, 309 (2002). 

9 Y. Alhassid and T. Rupp, Phys. Rev. Lett. 91, 056801 
(2003). 

10 I. S. Burmistrov, Y. Gefen, and M. N. Kiselev, JETP Lett. 
92, 179 (2010). 

11 I. S. Burmistrov, Y. Gefen, and M. N. Kiselev, Phys. Rev. 
B 85, 155311 (2012). 

12 R. L. Stratonovich, Dokl. Akad. Nauk SSSR 115, 1097 
(1957) [Sov. Phys. Dokl. 2, 416 (1958)]. 

13 J. Hubbard, Phys. Rev. Lett. 3, 77 (1959). 

14 J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 
108, 1175 (1957). 

15 J. von Delft and D. C. Ralph, Phys. Rep. 345, 61 (2001). 

16 P. W. Anderson, J. Phys. Chem. Solids 11, 26 (1959). 

17 A. Di Lorenzo, R. Fazio, F. W. J. Hekking, G. Falci, 
A. Mastellone, and G. Giaquinta, Phys. Rev. Lett. 84, 550 
(2000). 

18 G. Falci, R. Fazio, F. W. J. Hekking, and A. Mastellone, 
J. Low Temp. Phys. 118, 355 (2000). 

19 Y. Volokitin, J. Sinzig, L. J. de Jongh, G. Schmid, M. 
N. Vargaftik, and I. Moiseev, Nature (London) 384, 621 

(1996) . 

20 B. Muhlschlegel, D. J. Scalapino, and R. Denton, Phys. 
Rev. B 6, 1767 (1972). 

21 Y. Alhassid and J. Zingman, Phys. Rev. C 30, 684 (1984). 

22 B. Lauritzen, P. Arve, and G. F. Bertsch, Phys. Rev. Lett. 
61, 2835 (1988). 

23 A. K. Kerman and S. Levit, Phys. Rev. C 24, 1029 (1981). 

24 A. K. Kerman, S. Levit, and T. Troudet, Ann. Phys. (NY) 
148, 436 (1983). 

25 G. Puddu, P. Bortignon, and R. Broglia, Ann. Phys. (NY) 
206, 409 (1991). 

26 B. Lauritzen, A. Anselmino, P. Bortignon, and R. Broglia, 
Ann. Phys. (NY) 223, 216 (1993). 

27 R. Rossignoli and N. Canosa, Phys. Lett. B 394, 242 

(1997) . 



28 H. Attias and Y. Alhassid, Nucl. Phys. A 625, 565 (1997). 

29 A. L. Goodman, Nucl. Phys. A 352, 30 (1981). 

30 R. Rossignoli, N. Canosa, and P. Ring, Phys. Rev. Lett. 
80, 1853 (1998). 

31 R. Balian, H. Flocard, and M. Veneroni, Phys. Rep. 317, 
252 (1999). 

32 G. Falci, A. Fubini, and A. Mastellone, Phys. Rev. B 65, 
140507 (2002). 

33 K. Van Houcke, S. M. A. Rombouts, and L. Pollet, Phys. 
Rev. B 73, 132509 (2006). 

34 Y. Alhassid, L. Fang, and S. Schmidt, 
|arXiv:cond-mat/0702304| 

35 R. W. Richardson, Phys. Lett. 3, 277 (1963). 

36 R. W. Richardson, Phys. Rev. 159, 792 (1967). 

37 M. Schechter, Y. Imry, Y. Levinson, and J. von Delft, Phys. 
Rev. B 63, 214518 (2001). 

38 V. N. Gladilin, V. M. Fomin, and J. T. Devreese, Phys. 
Rev. B 70, 144506 (2004). 

39 Z. J. Ying, M. Cuoco, C. Noce, and H. Q. Zhou, Phys. 
Rev. B 74, 012503 (2006). 

40 S. Schmidt, Y. Alhassid, and K. Van Houcke, Europhys. 
Lett. 80, 47004 (2007). 

41 G. Falci, R. Fazio, and A. Mastellone, Phys. Rev. B 67, 
132501 (2003). 

42 K. Van Ho ucke, Y. A lhassid, S. Schmidt, and S. M. A. 
Rombouts, larXiv:1011.542T1 

43 M. Schechter, Phys. Rev. B 70, 024521 (2004). 

44 Y. Alhassid, S. Liu and H. Nakada, Phys. Rev. Lett. 99, 
162504 (2007). 

45 Y. Alhassid, H. A. Weidenmuller, and A. Wobst, Phys. 
Rev. B 72, 045318 (2005). 

46 S. D. Berger and B. I. Halperin, Phys. Rev. B 58, 5213 
(1998). 

47 If particle-number projection is carried out after the func- 
tional integral (|31[) is calculated, large errors are almost 
unavoidable because the canonical projection involves an 
oscillating sum or an integral of quantities whose only ap- 
proximate values are known. For the same reason, S z pro- 
jection is done for each value of Ao rather than after cal- 
culating the integral. 

48 Y. Alhassid, G. F. Bertsch, L. Fang, and S. Liu, Phys. Rev. 
C 72, 064326 (2005). 

49 P. Ribeiro, and A. M. Garcfa-Garci'a, Phys. Rev. Lett. 108, 
097004 (2012). 

50 I. S. Beloborodov, A. V. Lopatin, V. M. Vinokur, and K. 
B. Efetov, Rev. Mod. Phys. 79, 469 (2007). 



